function v=ppgrad(pp,xx)
%PPVAL  Evaluate piecewise polynomial derivative

if isstruct(xx) % we assume that ppval(xx,pp) was used
   temp = xx; xx = pp; pp = temp;
end

%  obtain the row vector xs equivalent to XX
sizexx = size(xx); lx = numel(xx); xs = reshape(xx,1,lx);
%  if XX is row vector, suppress its first dimension
if length(sizexx)==2&&sizexx(1)==1, sizexx(1) = []; end

% if necessary, sort xs 
if any(diff(xs)<0), [xs,ix] = sort(xs); end

% take apart PP
[b,c,l,k,dd]=unmkpp(pp);

% for each data point, compute its breakpoint interval
[ignored,index] = sort([b(1:l) xs]);
index = max([find(index>l)-(1:lx);ones(1,lx)]);

% now go to local coordinates ...
xs = xs-b(index);

d = prod(dd);
if d>1 % ... replicate xs and index in case PP is vector-valued ...
   xs = reshape(xs(ones(d,1),:),1,d*lx);
   index = d*index; temp = (-d:-1).';
   index = reshape(1+index(ones(d,1),:)+temp(:,ones(1,lx)), d*lx, 1 );
else
   if length(sizexx)>1, dd = []; else dd = 1; end
end

% ... and apply nested multiplication:
   v = c(index,1)*(k-1);
   for i=2:(k-1)
      v = xs(:).*v + c(index,i)*(k-i);
   end

v = reshape(v,d,lx);
if exist('ix'), v(:,ix) = v; end
v = reshape(v,[dd,sizexx]);
